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Abstract 


The  aim  of  this  project  was  to  develop  and  study  numerical  methods  for  the  solutions 
of  hyperbolic  equations  that  include  uncertainties.  In  particular  application  to  elec¬ 
tromagnetics  wa,ves  have  been  considered.  The  following  results  were  obtained  during 
this  research:  1.  The  nature  of  the  deterministic  systems  resulting  from  applying  the 
polynomial  chaos  method  to  an  hyperbolic  system  was  clarified,  nd  in  particular  the 
role  of  the  boundary  conditions.  2.  The  evolution  of  the  probability  density  func¬ 
tions  of  uncertainties  in  problems  described  by  hyperbolic  systems  was  found,  and  an 
application  to  the  Maxwells  equations  of  electromagnetics  had  been  carried  out.  3. 
Reduced  basis  methods  have  been  applied  to  several  wave  propagation  problems. 


1  Galerkin  Method  for  Wave  Equations  with  Un¬ 
certainty 


In  [10]  the  polynomial  chaos  method  applied  to  wave  problems  have  been  considered. 
In  recent  years  there  is  a  growing  interests  in  studying  efficient  numerical  methods  for 
solving  differential  equations  with  random  inputs.  Many  approaches  can  be  employed, 
among  which  the  Polynomial  Chaos  (PC)  based  methods  have  received  intensive  at¬ 
tention.  The  original  PC  method  was  developed  by  R.  Ghanern  and  was  inspired 
by  the  Wiener  chaos  expansion  which  uses  Herimte  polynomials  of  Gaussian  random 
variables  to  represent  random  processes.  Later  the  approach  was  extended  to  general¬ 
ized  Polynomial  Chaos  (gPC)  where  general  orthogonal  polynomials  are  adopted  for 
improved  representations  of  more  general  random  processes.  With  PC/gPC  serving 
as  a  complete  basis  to  represent  random  processes,  a  stochastic  Galerkin  projection 
can  be  used  to  transform  the  (stochastic)  governing  equations  to  a  set  of  determin¬ 
istic  equations  that  can  be  readily  discretized  via  standard  numerical  techniques. 
Although  such  a  Galerkin  approach  is  effective  in  many  problems,  its  application 
to  hyperbolic  problems  is  limited.  The  primary  reason  is  because  the  properties  of 
the  resulting  system  of  equations  from  a  Galerkin  projection  is  not  fully  understood. 
When  uncertainty  does  not  change  the  direction  of  the  characteristics,  the  Galerkin 
system  can  be  shown  to  be  hyperbolic  and  solved  in  a  straightforward  manner.  . 


We  studied  the  application  of  the  gPC  Galerkin  method  to  the  simulations  of  hy¬ 
perbolic  systems  that  contain  uncertainties.  In  general  these  uncertainties  may  enter 
through  initial  conditions,  boundary  conditions  or  through  uncertainties  in  the  co¬ 
efficients  of  the  problem.  In  solving  numerically  the  Maxwell’s  equations  of  electro¬ 
magnetics,  for  example,  one  often  faces  uncertain  material  properties  or  small  scales 
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that  can  not  be  resolved  and  have  to  be  modeled  statistically.  Here  we  deal  with  the 
case  that  the  coefficients  are  functions  of  random  variables.  In  particular  we  use  a 
scale  wave  equation  as  a  model  and  study  the  situation  in  which  the  inflow-outflow 
conditions  change  as  a  function  of  a  random  variable. 


We  have  shown  that  the  deterministic  system  is  a  symmetric  hyperbolic  system  with 
positive  as  well  negative  eigenvalues.  A  consistent  and  stable  method  of  imposing  the 
boundary  conditions  is  outlined.  The  boundary  conditions  are  not  satisfied  exactly 
at  the  boundaries  but  rather  to  the  order  of  the  scheme.  Convergence  of  the  scheme 
is  established. 


A  simple  scalar  equation  that  illustrates  the  difficulties  in  applying  the  (generalized) 
Polynomial  Chaos  to  hyperbolic  equations  is: 


du(x,  t,  y)  _  ,  .  d-u(x,t,y ) 

—  c{y)~ 


re  (-1,1),  t>  0, 


(1) 


at  dx 

where  c(y)  is  a  random  transport  velocity  of  a.  random  variable  j/Gfi  in  a  properly  de¬ 
fined  complete  random  space  with  event  space  Cl  and  probability  distribution  function 
p(y).  With  this  the  expectation  of  a  given  function  is  E[f(y))  —  f  f(y)p(y)dy. 


The  initial  condition  is  given  by 

u(x,Q,y)  =  u0(x,y).  (2) 


The  boundary  conditions  are  more  complicated  as  they  depend  on  the  sign  of  the 
random  transport  velocity  c(y).  A  well  posed  set  of  boundary  conditions  is  given  by: 

“(Miff)  =  up(t,y)  c(y)  >  0, 

u(-l,t,y)  =  uL(t,y)  c(y)  <  0. 

Equations  (  !)  --( 3)  complete  the  setup  of  the  problem. 


For  simplicity  we  will  discuss  the  case  of  random  variable  y  with  beta  distribution 
in  (-1,1)  (upon  proper  scaling).  In  this  case  the  expansion  functions  Pk  are  the 
(normalized)  Jacobi  polynomials.  (Note  this  includes  the  special  case  of  Legendre 
polynomials  with  uniformly  distributed  random  variable  y.) 

In  the  gPC  Galerkin  method  we  seek  an  approximation  to  the  true  solution  via  a 
finite-term  gPC  expansion 

N 

v(x,  t,y)  =  Y l  vk(x,  t)Pk(y )  (3) 

fc=o 
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and  project 


dvj^tyy)  _  dvjxXy)  _ 
dt  ^  dx 

onto  the  subspace  spanned  by  the  first  (JV  + 1)  gPC  basis  polynomials  and  obtain  the 
following  system 


dt 


dvk(x,t) 

dx 

k— 0 


j  =  0, . . . ,  JV. 


(4) 


If  we  denote  by  A  the  (N  +  1)  x  (N  +  1)  matrix  whose  entries  are  {cij,k}a<j,k<N  and 
v  =  (f’o,  •  •  •  ,  vn)t  a  vector  of  length  (N  +  1),  then  system  (4)  can  be  written  as 

dv(x,t)  _  dv(x,t) 

~t ~  =  k~ET‘  (5) 

Note  that  from  the  definition  ajk  =  c ikj ,  i.e.,  A  —  AT,  the  system  (5)  is  therefore 
symmetric  hyperbolic,  this  is  consistent  with  the  fact  that  the  original  equation  (1)  is 
hyperbolic  for  each  realization  of  y. 


A  less  trivial  question  is  the  nature  of  the  inflow-outflow  boundary  conditions.  The 
boundary  conditions  for  the  original  scalar  equation  (1)  depend  on  the  particular 
realization  of  the  random  variable  y  (see  (3)).  However  upon  the  Galerkin  projection 
in  the  random  dimension  the  deterministic  system  (5)  is  independent  of  y.  In  Theorem 
1  we  investigate  how  the  inflow-outflow  conditions  are  reflected  in  the  system  (5). 


Theorem  1: 


Consider  the  deterministic  system  (5)  where  the  coefficients  ajk are  defined.  Then  if 
c(y)  >  0  (reps.  c(y)  <  0)  for  all  y,  then  the  eigenvalues  of  A  are  all  non-negative 
(reps,  non-positive);  if  c(y)  changes  sign,  i.e.,  c(y)  >  0  for  some  y  and  c(y)  <  0  for 
some  other  y,  then  A  has  both  positive  and  negative  eigenvalues. 


Based  on  Theorem  1  have  developed  a  stable  way  of  employing  boundary  conditions 
for  the  GPC  method  for  hyperbolic  systems  of  equations.  A  detailed  discussion  is 
presented  in  [10] 
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2  Evolution  of  PDF  forHyperbolic  Equations  with 


Uncertainties 


In  [12]  the  evolution  of  the  probability  density  function  was  considered  for  PDEs  that 
simulate  wave  phenomena  containing  uncertainties  . 


Hyperbolic  equations  (such  that  the  Euler  equations  of  gas  dynamics  or  the  Maxwell’s 
equations  of  electromagnetics)  may  contain  uncertainties  in  the  coefficients  of  the 
equations  or  in  the  data  (initial  or  boundary  conditions).  One  way  of  treating  these 
uncertainties  is  to  model  them  by  random  variables  with  some  known  distributions. 
However,  even  if  the  initial  distributions  are  known,  the  probabilistic  properties  of  the 
solution  changes  in  time.  We  propose  to  study  this  behavior  in  order  to  understand 
the  time  development  of  the  solutions. 


In  [12]  the  cumulative  (or  probability)  distribution  function  (CDF)  of  the  output 
solutions  of  some  one-dimensional  wave  equations  fluctuating  by  the  random  changes 
in  the  media  parameters  was  computed.  We  will  review  here  this  work. 


The  numerical  tool  to  handle  randomness,  we  employ  the  technique  of  the  so-called 
polynomial  chaos  (PC)  expansions  The  main  advantage  of  the  PC  expansions  is  that 
we  separate  the  deterministic  and  the  non-deterministic  parts  in  the  functions  u  so 
that  we  can  compute  numerically  each  PC  mode  uk(x,t );  we  only  have  to  compute 
the  PC  modes  once  instead  of  evaluating  the  output  solution  for  each  £*. 


We  work  with  both  the  exact  solutions  and  the  numerical  ones  of  the  one-dimensional 
problems.  We  then  compare  the  statistics  for  both  solutions  (e.g.  PDF,  CDF,  mean, 
variance  and  moment  generating  function,  e.t.c.). 


Although  in  general  the  explicit  solutions  of  PDE  are  not  available  or  of  some  com¬ 
plicated  forms,  some  known  solutions  can  validate  the  numerical  (or  experimental) 
solutions.  Furthermore,  we  can  derive  explicitly  the  probability  distribution /density 
functions  of  output  solutions  and  we  can  then  trace  down  how  the  randomness  in  the 
input  data  affects  the  output.  We  also  derive  the  moment  generating  function,  mean 
and  variance. 


We  considered  the  1-D  Maxwell  equations  for  (u,v)T  =  (u^(xyt),v^(x,t))T  in  Q  = 
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(—oo,  oo)  x  (0,  oo): 


( 


\  V  /  t  \  € 

=  hi(x), 
v(x,Q)  =  /l2(®), 


u 

V 


in  fi, 


(2.1) 


where  (permeability,  permittivity)  =  (m,  e)  =  (/xj,  ei)  for  x  <  C  and  =  (/i2,  e2)  for  x  > 
C,  Mi>  M2>  £i,  t2  >  0,  C  is  a  random  variable  and  has  a  uniform  PDF  g(C)  =  ^X[-i,i](£)- 
Let  c  =  l/yfo/f,  Ci  =  l/i/eiMi  and  c2  =  l/1/e2/i2  with  Cj  /  c2. 


We  found  the  distribution  function  F(,s)  for  the  output  u.  By  definition,  we  have 

3 

F(s)  =  F(u  <  s)  =  P(u  <  s,  (x,  t )  €  Sd j ) .  (2.2) 

3=0 


We  can  deduce  the  distribution  function  F.  For  example,  if  hi(x)  =  x,  h2(x)  =  2x, 
Mi  —  £i  =  1)  M2  =  o2  =  1/3  (thus  Ci  =  1,  c2  =  3),  we  then  find  the  solutions: 


x  T  2t 

in  fio, 

4x  +  5 1 

-  3C 

in  f2i, 

^x  T  5t 

in  fi2, 

x  +  6 1 

in.  ^3. 

(2.3) 


At  the  point,  e.g.  (x,t)  =  (0, 1),  since  P( 2  <  s,  C  >  1)  =  P(6  <  s,  C  <  -3)  =  0,  we 
deduce: 


(2.4) 

(2.5) 

We  can  define  the  moment  generating  function  for  u,  m(t ),  and  apply  the  Monte 
Carlo  integration  too. 

oo  -i  M 

etuf(u)du~  —J2etu'-  (2.6) 

'°°  1  i=i 

Then  the  nth  moment  is: 

1  M 

ro(n>(0+)~  —  £(«*)"■  (2.7) 


F(s)  =  P(5  -  3C  <  s,  0  <  C  <  1)  +  P( 5  -  |  <  s,  -3  <  C  <  0) 


i 

V  2 


if  s  <  2 
if  2  <  s  <  5 
if  s  >  5 


+ 


0 

3(5  —  5) 
i 2 

2 


if  s  <  5 
if  5  <  s  <  f 

if  S  >  y 


5 


The  mean  and  variance  can  be  evaluated  as  follows. 


E(u )  =  m/(0+)  ~ 


i—  1 


Var(u)  =  m"( 0+)  -  ( E(u ))2  - 


1 

M 


M 


£M2  - 


(2.8a) 


(2.8b) 


3  Certified  Reduced  Basis  Methods 


In  [1]  we  proposed  certified  reduced  basis  methods  to  enable  the  efficient  and  reliable 
evaluation  of  a  general  output  that  is  implicitly  connected  to  a  given  parameterized 
input  through  the  harmonic  Maxwell’s  equations.  The  truth  approximation  and  the 
development  of  the  reduced  basis  through  a  greedy  approach  is  based  on  discontinuous 
Galerkin  approximations  of  the  linear  partial  differential  equation.  The  development 
allows  the  use  of  different  approximation  spaces  for  solving  the  primal  and  the  dual 
truth  approximation  problems  to  respect  the  characteristics  of  both  problem  types, 
leading  to  an  overall  reduction  in  the  off-line  computational  effort. 


The  main  features  of  the  method  are:  i)  rapid  convergence  on  the  entire  set  of  pa¬ 
rameters,  ii)  rigorous  a  posteriori  error  estimators  for  theoutput  and  iii)  a  parameter 
independent  off-line  phase  and  a  omputationally  very  efficient  on-line  strategy  to  en¬ 
able  the  rapid  solution  of  many  query  problems  arising  in  control,  optimization,  and 
design.  The  versatility  and  the  performance  of  this  approach  is  shown  through  some 
numerical  experiments,  illustrating  the  modeling  of  material  variations,  problems  with 
resonant  behavior,  and  applications  related  to  radar  scattering  prediction. 


In  [2]  A  Galerkin  approach  was  combined  with  a  reduced  basis  method  for  the  eval¬ 
uation  of  outputs  of  interest  implicitly  depending  on  a  given  input  via  the  resolution 
of  a  PDE  issue  from  a  harmonic  wave  propagation  problem.  The  main  features  of 
the  method  are:  i)  rapid  convergence  on  the  entire  set  of  parameters,  ii)  a  posteriori 
error  estimators  for  the  output  and  iii)  an  off-line  (parameter  independent)  on-line 
(very  fast)  computational  strategy.  In  the  present  paper  we  allow  the  use  of  different 
approximation  spaces  for  solving  the  primal  and  dual  truth  approximation  problems 
reducing  the  off-line  computational  effort. 

In  [3]  In  a  posteriori  error  analysis  of  reduced  basis  approximations  to  affinely  pa¬ 
rameterized  partial  differential  equations,  the  construction  of  the  lower  bounds  for 
coercivity  and  inf-sup  stability  constants  is  required.  In  [1]  the  authors  presented  an 
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efficient  method  compatible  with  an  off-line/on-line  strategy,  where  the  online  compu¬ 
tation  is  reduced  to  minimizing  a  linear  functional  (assuming  the  affine  dependence 
on  the  parameters)  under  a  few  linear  constraints  reflecting  stability  information. 
These  constraints  depend  on  nested  sets  of  parameters  obtained  iteratively  using  a 
greedy  algorithm.  In  the  paper,  we  improved  this  method  so  that  it  becomes  more 
efficient  and  have  a  nice  property,  namely,  the  computed  lower  bound  is  monotoni- 
cally  increasing  with  respect  to  the  nested  sets.  The  new  method  is  applied  to  the 
construction  of  the  inf-sub  stability  constant  of  an  electromagnetic  cavity  problem, 
providing  good  lower  bounds  and  capturing  the  resonance  lines  very  efficiently  and 
accurately. 


4  Problems  with  random  inputs 


In  [5]  we  proposed  a  way  of  accounting  for  the  lack  of  detailed  knowledge  about 
material  shapes  in  computational  time-domain  electromagnetics.  We  use  Legendre- 
Gauss-Lobatto,  Stroud-2  and  Stroud-3  quadrature  formulae  to  solve  the  resulting 
stochastic  equation  and  we  show  the  efficiency  of  the  proposed  method  over  statistical 
Monte-Carlo  simulations.  We  also  show  how  the  radar-cross-section  in  scattering  is 
affected  by  the  uncertainty  in  shape  of  the  objects  and  by  the  direction  of  the  incident 
field. 


In  CHL  we  discuss  computationally  efficient  ways  of  accounting  for  the  impact  of 
uncertainty,  e.g.,  lack  of  detailed  knowledge  about  sources,  materials,  shapes,  etc.,  in 
computational  time-domain  electromagnetics.  In  contrast  to  classic  statistical  Monte 
Carlo  based  methods,  we  explore  a  probabilistic  approach  based  on  high-order  accu¬ 
rate  expansions  of  general  stochastic  processes.  We  show  this  to  be  highly  efficient 
and  accurate  on  both  one-  and  two-dimensional  examples,  enabling  the  computa¬ 
tion  of  global  sensitivities  of  measures  of  interest,  e.g.,  radar-cross-sections  (RCS)  in 
scattering  applications,  for  a  variety  of  types  of  uncertainties. 


In  [7]  Recently  there  has  been  a  growing  interest  in  designing  efficient  methods  for 
solutions  of  ordinary/partial  differential  equations  with  random  input.  The  thisend, 
Stochastic  Galerkin  (SC)  methods  appear  to  be  superior  to  other  non-sampling  meth¬ 
ods,  and  in  many  cases,  to  several  sampling  methods.  However,  when  the  governing 
equations  take  complicated  forms,  numerical  implementations  of  SG  methods  can 
become  non-trivial  and  care  is  needed  to  design  robust  and  efficient  solvers  for  the  re¬ 
sulting  equations.  On  the  other  hand,  the  traditional  sampling  methods,  e.g.,  Monte 
Carlo  methods,  are  straightforward  to  implement,  but  they  do  not  offer  convergence 
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as  fast  as  stochastic  Galerkin  methods.  In  this  paper,  a  high-order  Stochastic  collo¬ 
cation  approach  is  proposed.  Similar  to  SG  methods,  the  collocation  methods  also 
take  advantage  of  the  smoothness  assumption  of  the  solutions  in  random  space  to 
achieve  fast  convergence.  One  the  other  hand,  the  numerical  implementation  of  the 
stochastic  collocation  is  trivial  as  it  only  requires  repetitive  runs  of  an  existing  de¬ 
terministic  solver,  similar  to  Monte  Carlo  methods.  The  computational  cost  of  the 
collocation  methods  depends  on  the  choice  of  the  collocation  points  and  we  present 
several  feasible  constructions.  One  particular  choice,  based  on  sparse  grids,  depends 
weakly  on  the  dimensionality  of  the  random  space  and  is  more  siutable  for  high  ac¬ 
curacy  computations  of  practical  applications  with  large  dimensional  random  input. 
Numerical  examples  are  presented  to  demonstrate  the  accuracy  and  efficiency  of  the 
stochastic  collocation  methods. 


In  [8]  It  is  well  known  that  the  steady  state  of  an  isentropic  flow  in  a  dual-throat 
nozzle  with  equal  throat  areas  is  not  unique.  In  particular  there  is  a  possibility 
tha  tthe  flow  contains  a  shock  wave,  whose  location  is  determined  solely  by  the 
initial  condition.  In  this  paper,  we  consider  cases  with  uncertainty  in  this  initial 
condition  and  use  generalized  polynomial  chaos  methods  to  study  the  steady-state 
solutions  for  stochastic  initial  conditions.  Special  interest  is  given  to  the  statistics 
of  the  shock  location.  The  polynomial  chaos  (PC)  expansion  modes  are  shown  to 
be  smooth  functions  of  the  spatial  variable  x,  although  each  solution  realization  is 
discontinuous  in  the  spatial  variable  x.  When  the  variance  of  the  initial  condition  is 
small,  the  probability  density  function  (PDF)  of  the  shock  location  is  computed  with 
high  accuracy.  Otherwise,  many  terms  are  needed  in  the  PC  expansion  to  produce 
reasonable  results  due  to  the  slow  convergence  of  the  PC  expansion,  caused  by  non¬ 
smoothness  in  random  space. 
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